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Abstract. The modified gravity is considered to be one of possible explanations of the 
accelerated expansions of the present and the early universe. We study effects of the modified 
gravity on big bang nucleosynthesis (BBN). If effects of the modified gravity are significant 
during the BBN epoch, they should be observed as changes of primordial light element abun¬ 
dances. We assume a f{G ) term with the Gauss-Bonnet term G , during the BBN epoch. A 
power-law relation of df /dG oc t p where t is the cosmic time was assumed for the function 
f(G ) as an example case. We solve time evolutions of physical variables during BBN in the 
f(G ) gravity model numerically, and analyzed calculated results. It is found that a proper 
solution for the cosmic expansion rate can be lost in some parameter region. In addition, 
we show that calculated results of primordial light element abundances can be significantly 
different from observational data. Especially, observational limits on primordial D abundance 
leads to the strongest constraint on the f(G) gravity. We then derive constraints on param¬ 
eters of the f(G) gravity taking into account the existence of the solution of expansion rate 
and final light element abundances. 
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1 Introduction 

The accelerated expansion of the present universe has been verified by observational data on 
magnitude-redshift of type la supernovae [1, 2], In standard cosmology, this acceleration is 
described by the cosmological constant or the A term in the gravitational action. Dark energy 
or modified gravity theory are also widely considered to explain the late time acceleration and 
one of the challenges is to discriminate the dark energy model and the modified gravity theory 
through the observations. For example, each of those models shows the different histories of 
cosmic expansion (e.g., [3]) and growth rate for the large scale structure (e.g., [4, 5]). 

One of feasible models of the modified gravity includes the Gauss-Bonnet (GB) correc¬ 
tion, G = R 2 — AR^R^ + R pupa R ,1L,pa , which is a simple extension of Einstein gravity. The 
GB term does not produce any ghost particles as well as any problems with unitarity [6]. 
Additionally, the equations of motion do not contain higher than second order in temporal 
derivatives, so there does not exist any instability problem. We generalize the GB correction 
as /(G) in which /(G) is a general function of G. Cosmological effects of this /(G) model 
have been investigated intensively as a simple extension of the general relativity. In general, 
modifications of the general relativity are restricted from observed celestial motions in the 
solar system. However, no correction to the Newton’s law and no instability are induced at 
the present universe in the /(G) models when the functional shape of /(G) is fine-tuned [3]. 
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Therefore, the f(G) gravity can escape from the constraint from the Newton’s law, so that it 
can be a candidate theory for the accelerated expansion of the universe. 

It has been suggested that epochs realized in the standard ACDM model can be described 
by some functions f{G ) in the f(G) models and that additional degrees of freedom in the 
f(G) model can be tested in future observations of cosmic expansion (e.g., [7]). Constrains on 
modified gravity models are, however, derived not only from the cosmic expansion rate but 
also from considerations of the cosmological perturbation (e.g., [4, 5]). The latter constraint 
has been found to be very strong. It has been argued that it is theoretically interesting to 
check whether viable f(G) models satisfy the weak energy condition, i.e., T^ lu U tJ 'U L ' > 0 for 
all timelike vectors U p with T /1U the stress-energy tensor [8]. 

The modified gravity can also lead to a change in cosmic expansion rate in the early 
universe. It, therefore, provides a solution to the horizon problem, flatness problem, and 
other problems related to observations of cosmic microwave background radiation [9]. In this 
way, the modified gravity is one of mechanisms for changing the cosmic expansion rate in the 
very early and the present epochs. It is, however, possible that effects of the modified gravity 
will be detected in cosmological observables which look to be consistent with astronomical 
observations at the moment. If the cosmic expansion rate is different from that in the standard 
cosmological model during the big bang nucleosynthesis (BBN) epoch, primordial abundances 
of light elements can be different from those predicted in standard BBN (SBBN) model. These 
differences may be detected in future observations of elemental abundances. Therefore, it is 
meaningful to investigate effects of the modified gravity on BBN theoretically, and predict 
light element abundances in the modified gravity models. We can also derive constraints on 
the modified gravity during the BBN epoch independently of those during the extremely early 
and present epoch. Effects of the f(R) gravity, where f(R) is a general function of R, i.e., 
Ricci scalar, has been studied for a specific type of f(R) oc R n as a simple extension of the 
general relativity [10-12]. Primordial abundances of D, 3,4 He, and 6, 'Li have been calculated 
in the model, and a constraint on the model has been derived from comparison of calculated 
abundances with observational data [12]. 

In this paper, we analyze effects of the f(G ) modified gravity during BBN epoch, and 
derive constraints on the model based on observed light element abundances. In section 2, 
the f{G ) model is introduced. We study a specific case of £(G) = df(G)/dG oc t p with t 
the cosmic time and p a power-law index since numerical BBN calculations in the model are 
easily performed. In section 3, our BBN code and treatment of modified gravity effects are 
explained. In section 4, observational constraints on primordial light element abundances 
are described. In section 5, we show results of BBN in the f(G) gravity model, and derive 
constraints on the model. In section 6, we summarize this study. In this paper, we adopt 
natural units of c = fee = 1, where c is the light speed and k& is the Boltzmann constant, 
and the convention of [V Q , Vg\A IJ ' = R IJ ‘ l/a pA 1/ and R^ u = R x ^\ v . 

2 Model 

2.1 f(G) gravity 

Equations of motion in the f(G) gravitational model are shown in this section. The action 
taken in this paper is given by 

S= 2 K 1 / [ i? + k2 ^( G )] + S m(g!J,wAm), (2-1) 
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where k 2 = 87 tQ is defined with Q the Newton’s constant, g iW is the metric tensor, g is the 
determinant of the metric tensor, S m is the action of the matter field <p m . We assume the 
spatially flat Friedmann-Lemaitre-Robertson-Walker metric as the same as in the standard 
cosmological model 

ds 2 = — dt 2 + a(t) 2 ( dx 2 + dy 2 + dz 2 ) , ( 2 . 2 ) 

where a(t ) the scale factor of the universe. For this metric, the nonzero components of the 
Ricci tensor are given by 


The Ricci scalar is given by 


The GB term is given by 


Roo — “3-, 

a 


Rij = a 


+ 2 


Jiy 


R = 6 



G = R 2 - 4R I , IJ R> W + R^R^ 
= 24 H 2 + FT 2 ) . 


(2.3) 

(2.4) 


(2.5) 


( 2 . 6 ) 


For a matter, on the other hand, we assume a perfect fluid described with an energy 
density p(t) and pressure p(t) 

= diag (-p, p,p,p). (2.7) 

The held equation for the f{G ) gravity is then derived by varying the action (eq. (2.1)) with 
respect to the metric tensor, 

- 9 -fR + R vv + {-^ff(G) + 2 (R lw + D IW ) ( f'R) 

+8 V a [V„ (f'Rvx) + V, (f'Rpx)] 

if Rap) ~ 4D (f'Rpu) + 2 f'Rp a R a v 
+2 f'Rf^RvaPi + 4V p V cr (f'R fipi/a 4“ f Rjimjp) J — K Tjiv 1 ( 2 - 8 ) 


where f'(G ) = df(G)/dG is the derivative with respect to G , is the covariant derivative 

operator, □ = p /i '"V^V,, = V P V M is the D’Alambertian operator, and the differential tensorial 
operator Dp U is defined as 

Dp U = pS7 u . (2.9) 

Especially, varying the action with respect to the 0-0 and i-i components of metric tensor, 
one derives 


f(G) - Gf 


+ 12 H 3 f 


f(G) - Gf , , „ 2 


+ 4 H z f + 2 Hf + 8 HHf 


3 H 2 = k 2 P 

( 2 . 10 ) 

3 H 2 = -n 2 p, 

( 2 . 11 ) 


where H = a/a is the expansion rate of the universe. 
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In addition, the conservation of the stress-energy tensor V ,L T fw holds, which results in 
the equation 

p + 3H(p + p) = 0. (2.12) 

This equation can be found using eqs. (2.10) and (2.11) also. 


2.2 Assumption on f{G) 

In this paper, we constrain the model space of f{G). We define a variable £ = f, and derive 
equations from eqs. (2.10) and (2.11), 


H(H-,£,p,p) 


2 k 2 P - 6 H 2 - 24 H 3 K 2 i - K 2 f{G) + k 2 £G = 0 


ft 2 

4 H 2 


- 

Ip + p) 

2 

\ + 4 K 2 Hi 



We assume that the £ term scales as a power-law function of time, i.e., 


(2.13) 

(2.14) 


£(*) = £ 0 (t/t 0 ) p , 


(2.15) 


where £(i) is dimensionless, £o = £(to) is dimensionless parameter for the £ value at to = 1 
s, and p is a power-law index. This model is a toy model for time-dependent f(G) term. 
In this paper, we show how to solve BBN in the f{G ) modified gravity model exactly in 
numerical calculation. Usually, detailed BBN calculations are performed with time taken as 
the evolution parameter. Because of this assumption of the explicit time dependence of £, 
values of £, £, and £ can be specified before solving H(t), H(t), G(t ), and f(G) using eqs. 
(2.6), (2.13), (2.14), and (3.2). Since we do not need to treat the parameter £ as unknown 
in solving the system of equations, it becomes somewhat simple to solve physical variables in 
the present model. It would be possible to solve BBN in other models of f(G) gravity based 
on the method described below. However, methods of the solution can be more complicated 
depending on the adopted models, and it is beyond the scope of this study to investigate 
many possible models of f(G ) function. 

Effects of the £ term on primordial light element abundances are constrained from 
observations of the primordial abundances. Since the primordial abundances are sensitive 
to the cosmic expansion history only in the BBN epoch, the GB term during BBN can 
be constrained. We then assume that the GB term exists in the temperature range of 
100 >T 9 = T/(10 9 K) > 0.01. 


3 BBN calculation with modified expansion rate 

The public BBN calculation code [13, 14] is utilized and modified. In this study, the ef¬ 
fective number of neutrino species is assumed to be three. We updated reaction rates of 
nuclei with mass numbers < 10 using the JINA REACLIB Database [15] (the latest version 
taken in December, 2014). The neutron lifetime is the central value of the Particle Data 
Group, 880.3 ±1.1 s [16]. The baryon-to-photon ratio is taken from the accurate value, 
(6.037 ± 0.077) X 10~ 10 [17], corresponding to the baryon density in the base ACDM model 
(Planck+WP±highL±BAO) determined from Planck observation of cosmic microwave back¬ 
ground, Q m h 2 = 0.02205 ± 0.00028 [18]. 

In general, two independent equations among equations of motion are used in BBN 
numerical codes. For example, in the Kawano’s code, equations of the Hubble expansion rate 
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H and the time derivative of temperature dT/dt are used. In the present modified gravity 
model, the Hubble expansion rate is given by eq. (2.13). The energy conservation equation, 
i.e., eq. (2.12), is, however, the same as that in the SBBN model. We can, therefore, use the 
same equation for the time evolution of temperature as that in the SBBN model (eq. (D.26) 
in ref. [13]). Then, only one modification of the Hubble rate should be added to the code for 
BBN network calculations. We should solve the Hubble expansion rate H(t) and H(t), and 
calculate the function 

f(G) = f £(t)dG. (3.1) 

JGi 

We assume that the £ term exists in the temperature range of Tg = [10 2 ,10 -2 ]. This 
temperature range corresponds to the time range of t = 0( 10 ~ 2 — 10 6 ) s in SBBN model 
although the time range depends on the model parameters in the f(G) model. The most 
important temperature relevant to SBBN is T 9 ~ 1. The power-law of £(i) = £o(t/to) p is 
assumed with the range of index p = [—2,6]. This range is chosen since it includes critical 
values of p = 2 and 4, and effects of f(G ) gravity outside the range are trivial (see section 
5). Both of the positive and negative £ cases are considered. Another parameter is the initial 
value of the f(G) term, i.e., f(Gj). The initial time corresponding to the initial temperature 
of Tg = 10 2 is assumed to be the same as that of SBBN. This choice of initial time only affects 
the normalization of the £ value and is therefore not important. 

We adopt the Newton-Raphson method to search for a correct solution with careful 
attention to the possibility of finding fake solutions. We used the following technique. For a 
given time, the temperature is calculated, and various physical quantities including the the 
energy density and pressure are derived. Since the function £(f) has been given (eq. (2.15)), 
the H value is given as a function of H (eq. (2.14)). The GB term G is then given as a 
function of H also (eq. (2.6)). The function f(G ) is calculated by integration of £(t) with 
respect to the GB term (eq. (3.1)). The value of G at time t + At is estimated with the 
equation, 

f(G(t + At)) = f(G(t)) + ^-- [G(t + At) — G(t)], (3.2) 

where f(G{t)) and £(f) is the values of f(G) and £ at the previous time in the calculation. 
The function g(H;^,p) is then given as a function of H , and a solution of g(H;£,p) = 0 is 
searched for. In the numerical calculation, the following equations are utilized. 


dg 

dH 

dG 

dH 

dH 

dH 

df(G) 

dH 


-12 H + k 2 


—72if 2 £ + £^777 - 


dG df(G) 


96 H 3 + 48 HH + 24 H 2 ^- 


2k 2 3FT 2 £ - 2H , 




dH 
d 
dH 


dH 


4 k 2 £ 


1 + 4 K 2 iT£ 1 + 4 K 2 iT£ 

£(t) + £(t + At) dG 
2 ~dH 


-H 


(3.3) 

(3.4) 

(3.5) 

(3.6) 


At the initial time in the calculation, the value for an initial guess of Hubble rate is 
given by H gVLess = HsbbNj where F^sbbn is the rate in SBBN model at the same temperature. 
At later times, the initial guesses are given by solutions found at the previous times in the 
calculation. The equation of the Hubble rate (eq. (2.13) with eqs. (2.6), (2.14), and (3.2)) 
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is a transformed quintic equation. Their coefficients, however, evolve as a function of time. 
Therefore, the number of real roots can change with time. It is then possible that a real 
root for the Hubble rate evolves to a complex root at some time. This root is of course 
inappropriate solution for the expansion rate of the real universe. Also if a root becomes 
negative at some time, it is inappropriate. Therefore, we need to evaluate whether there is a 
real positive solution for a parameter set of (£o,p, f(Gi)). We classify situations of the root 
finding into six cases: (1) a positive real root is found and it is not different from the previous 
root by more than three orders of magnitude, (2) any real root is not found, (3) a negative real 
root is found, (4) a positive real root is found but is is larger than that in SBBN by a factor 
of more than 10 s , (5) a positive real root is found but it is different from the previous root 
by more than three orders of magnitude, (6) a positive real root is found which is different 
from case (5), but the sign of df /dH is opposite to that at the previous time. 

We regard that an appropriate root is obtained only in the case (1). In cases (2)-(6), the 
solution most probably becomes negative, or complex, or too large so that they are included in 
a “no solution” region. The case (2) can be caused by the disappearance of the real solution. 
In this case, we do not have a smooth positive solution which continuously exists during 
the temperature range for the calculation. The case (3) can be caused by a decrease of the 
root value below zero. The case (4) can be caused by an unrealistically large root, and was 
excluded from real solutions. The case (5) can be caused by finding a root different from that 
in the previous time, or a change of the root value across zero. The case (6) can be caused 
by finding a different root. It indicates that one positive real root becomes negative, and the 
code find another positive real root next to it, most probably. 


4 Observed light element abundances 

Calculated BBN results are compared to the following observational constraints on light 
element abundances. 

The primordial 4 He abundance is estimated with observations of metal-poor extragalac- 
tic H II regions. We use the latest determination of Y p = 0.2551 ± 0.0022 [19]. When the 
central values of adopted reaction rates, the neutron lifetime, and the baryon-to-photon ratio 
are used, the calculated abundances in the SBBN model is out of the 2a observational limit. 
The 4a range is then adopted in this study. 

The primordial D abundance is estimated with observations of metal-poor Lyman-a 
absorption systems in the foreground of quasi-stellar objects. We use the weighted mean value 
of D/H= (2.53 ± 0.04) x 10 -5 [20], and adopt its 4er range since the 2a range is inconsistent 
with the theoretical abundances in the SBBN model, similarly to the 4 He abundance. 

3 He abundances are measured in Galactic H II regions through the 8.665 GHz hyperfine 
transition of 3 He + ion. These are not the primordial abundance but present values which 
have contributions from Galactic chemical evolution taking into account production and de¬ 
struction of nuclei in stars. Nevertheless, it is very hard to reduce elemental abundances 
significantly in standard Galactic chemical evolution theory. We then adopt the 2a upper 
limit from the abundance 3 He/H=(1.9 ± 0.6) X 10~ 5 [21] in Galactic H II regions, as a rough 
guide. 

The primordial 7 Li abundance is estimated with observations of Galactic metal-poor 
stars. We use the abundance log( 7 Li/H)= —12 + (2.199 ± 0.086) derived in a 3D nonlocal 
thermal equilibrium model [22], and adopt its 2 a range. 
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power law index p 


Figure 1. Boundary of regions with and without proper cosmological solutions in the parameter 
plane of (p, £o) in the case of £o > 0. For this figure the initial value of /(G) is fixed to be /(Gi) = 0. 
The region marked with ‘no solution’ above the boundary corresponds to parameter sets which have 
no proper solution of cosmic expansion. This region is therefore excluded. 


6 Li abundances in Galactic metal-poor stars have also been measured. We adopted the 
least stringent 2 a upper limit of all limits for stars reported in [23], i.e., 6 Li/H=(0.9± 4.3) x 
10~ 12 for the G64-12 (nonlocal thermal equilibrium model with 5 free parameters). 


5 Result 


5.1 Case of £o > 0 and /(Gi) = 0 

Figure 1 shows boundary of regions with and without proper cosmological solutions in the 
parameter plane of (p, £q) hi the case of £o > 0. The initial value of /(G) is fixed as 
/(Gi) = 0 for purposes of illustration. The region marked with ‘no solution’ above the 
boundary corresponds to parameter sets which lead to no proper solution of cosmic expansion 
in the BBN calculation. This region is therefore excluded. We find that the constraint in the 
parameter plane of (p, £q) f° r £o > 0 is predominantly given by the continuous existence of 
the solution H > 0 during the BBN epoch. Final abundances of light elements are within 
the adopted observational limits, in the parameter region below the boundary. The limits on 
light element abundances are, therefore, less constraining than the existence of a solution. 

The reason for the shape of the boundary is explained as follows. We consider deviations 
of the cosmic expansion rate from that in the standard cosmological model. First, some 
characteristics of the standard model are reviewed. The expansion rate during the radiation 
dominated epoch is given by H = a/a = 1/(2 1). Also the following equations holds, 


H 2 = 
H = 
H + H 2 = 


K 

(5.1) 

-y (P + P) 

(5.2) 

-yG* + 3p) 

0 

(5.3) 
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Since H + H 2 < 0 is always satisfied, the inequality G < 0 holds (eq. (2.6)). The time 
derivative of G is given by 

dG 4 a d 

s = - rdd pip+m - (5 - 4) 

Therefore, the inequality dG/dt > 0 is satisfied. 

If the deviation of expansion rate is not large, it follows that 

- 6 H 2 Z ~t~ 2 (5.5) 

-24 H 3 k 2 £ Z -H 3 t oc - 1~ 3 (pCo* p_1 ) = -pCo * p ~ 4 (5.6) 

\j-p- 41 1 

-K 2 f(G) Z - j idG (X - J (Hot?) d{-t~ A ) (X -Co J Hot p ~ 5 dt <x (5.7) 

k 2 £G ZHG oc (6 t p ) (-t~ A ) = -Co t p ~ A . (5.8) 

We note that once the expansion rate deviate from the value in the standard model signifi¬ 

cantly, the relations of H oc t ” 1 and G oc — t~ A are broken and the above scalings no longer 
hold. 

p < 4: As t increases, amplitudes of the second, third, and fifth terms after the first equality 
in eq. (2.13) decrease. However, the amplitude of the fourth term does not change much since 
the early time of t ~ t\ contributes to the integral (eq. (5.7)) predominantly. As a result, the 
fourth term becomes more and more important relatively. Because of f(G) > 0 for Co > 0 
(eq. (5.7)), the fourth term is negative and decelerates the cosmic expansion. The f(G) 
term, therefore, works similarly to a negative dark energy A. When this term becomes large, 
we lose a real positive root of the expansion rate which smoothly connects to the root in the 
standard model. In the parameter region marked as ‘no solution’, the solution is lost until the 
cosmic temperature decreases down to the final value of Tg = 10 -2 . On the boundary of the 
‘no solution’ region, solutions disappear right before the temperature decreases to Tg = 10~ 2 . 
On the boundary, therefore, the equation is satisfied, 

f(G) « 2p, (5.9) 

at Tg = 1CP 2 . This equation holds since the second, third, and fifth terms in eq. (2.13) are 
negligible for H —>• 0. When the power-law index p is larger, the amplitude of the fourth 
term oc | [t p_4 ] t . | ~ tf 4 is smaller. Therefore, the effect of the f(G) term is smaller and 
larger amplitudes of Co are allowed from observational constraints. This is the reason for the 
upward-sloping curve for p < 4. 

p > 4: As t increases, amplitudes of the third to fifth terms in eq. (2.13) increase, while 
that of the second term decreases. Since all of the third to fifth terms are negative, their 
sum works similarly to a negative A. The boundary then corresponds to the point of eq. 
(5.9) also for p > 4. All of these terms roughly scale as t p ~ A for f > f; (eqs. (5.6)-(5.8)). 
When the power-law index p is larger, the quantity t p ~ 4 is larger so that the effects of the 
three terms are larger. As a result, smaller amplitudes of Co are allowed. This results in the 
downward-sloping curve for p > 4. 

We note that a magnitude of the parameter Co allowed from the BBN constraint is 
typically huge. For example, we consider a case of p = 4 in which the cosmic expansion 
rate is not so much different from that in the standard cosmological model. The third to 
fifth terms in eq. (2.13) then become constants (see eqs. (5.6)-(5.8)). Since the three terms 



scale similarly, we just take the fifth term for simplicity. Then, the amplitude of £ which 
significantly affects the expansion rate is estimated as follows. In this case, the first, the 
second and the fifth terms in eq. (2.13) are of the same order of magnitude, i.e., 

2n 2 p ~ 6H 2 ~ k 2 |£||G|. (5.10) 


The GB term in the standard model is given by 


G = 2AH 2 [H + H 2 ^j ~ —|/cV, (5.11) 

where eqs. (2.6), (5.1), and (5.3) were used, and the radiation dominated epoch is assumed. 
The energy density in the radiation dominated universe is given by 

p=^g*T\ (5.12) 

where g* is the relativistic degrees of freedom for the energy density [24]. Inserting eqs. (5.11) 
and (5.12) into |£| ~ 2p/\G\ from eq. (5.10), we obtain 


l£l ~ 


3 1 

4 K A p 


45 

128vr 4 5 * 



(5.13) 


where mpi = Q~ 1 ^ 2 = \/&k/k = 1.22 x 10 19 GeV is the Planck mass. Since the factor 
(mpi/T) 4 is huge during the BBN epoch of T ~ 1 MeV, constraints on the present modified 
gravity model exclude huge amplitudes of £. 


5.2 Case of £o < 0 and /(G;) = 0 

5.2.1 Parameter search 

Figure 2 shows the same boundary (black line) as in figure 1 but in the parameter plane of 
(p, |£o|) for £o < 0. Contours of calculated light element abundances are also drawn. Solid 
and dashed lines for D (green line), 3 He (purple lines), and 4 He (red lines) correspond to the 
observational upper (‘high’) and lower (‘low’) limits, respectively, on their abundances. Blue 
solid and dashed lines marked with ’ 7 Li ‘obs’ correspond to observational upper and lower 
limits, respectively, on 'Li/H abundances. We find that observational limits on primordial 
abundances give important constraints in the region of p < 4, contrary to the case of £o > 0. 

The shapes of the solution boundary and the abundance contours are explained as fol¬ 
lows. Deviations of the cosmic expansion rate from that in the standard cosmology is consid¬ 
ered again. 

p<2: As t increases, the fourth term in eq. (2.13) gradually becomes dominant. Because 
of /(G) < 0 for £o < 0, the fourth term is positive (eq. (5.7)). Therefore, the term causes 
an acceleration of the cosmic expansion similarly to a positive A. When this term becomes 
large, the expansion rate deviates from the standard rate and the universe inflates. Since the 
fourth term scales as ~ — 1\~ , a larger power-law index p gives smaller amplitudes of the 
fourth term. The effect of the /(G) term is, therefore, smaller and larger |£q| are allowed. 
This fact explains the upward-sloping curves for the solution boundary and the abundance 
constraints for p < 2. For a fixed power-law index p, larger |£o| values causes larger H values 
or faster cosmic expansion. Effects of the /(G) gravity on light element abundances are, 
therefore, larger. Reasons of the complicated shape of the boundary are described below in 
section 5.2.3. 
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power law index p 


Figure 2. Same boundary (black line) as in figure 1, but in the parameter plane of (p, |£o|) for £o < 0. 
Contours of calculated light element abundances are also shown. Solid and dashed lines for D (green 
line), 3 He (purple lines), and 4 He (red lines) correspond to the observational upper (‘high’) and lower 
(‘low’) limits, respectively on their abundances. Blue solid and dashed lines marked with i7 Li obs’ 
correspond to observational upper and lower limits, respectively, on 'Li/H abundances. Four orange 
stars located in the ‘no solution’ region correspond to the parameter set analyzed in section 5.2.2. 


P <4: As t increases, amplitudes of components which are proportional to t p 4 in the 
third to fifth terms in eq. (2.13) become larger relative to the term —6 H 2 oc t ~ 2 . In a late 
time of BBN, then a large deviation of the expansion rate is induced, and a proper solution 
possibly disappears (see section 5.2.2). A larger power-law index p gives larger amplitudes 
of the third to fifth terms. Since the effect of modified gravity is thus larger, smaller |£q| 
values are allowed. For this reason, the solution boundary is downward-sloping in the region 
2 < p < 4. 

Jr r^j 

P > 4: As t increases, amplitudes of the third to fifth terms in eq. (2.13) become large 
relative to that of the second term. The third to fifth terms are positive so that their sum 
works similarly to a positive A. When the amplitude |£q| is larger than a critical value, a 
real positive root of H is lost (see section 5.2.2). All of these terms roughly scale as t p ~ 4 for 
t t[. When the power-law index p is larger, the quantity t p ~ A is larger and the effects of 
the three terms are larger. Smaller |£q| values are then allowed. This is the reason for the 
downward-sloping curve for p > 4, similarly to the case of £o > 0. 

5.2.2 Examples of cosmic evolution 

Evolutions of the expansion rate are illustrated for four parameter sets in the ‘no solution’ 
region, which are indicated by orange stars in figure 2. Since the boundary of the existence 
of proper solution has a complicated shape, we check the cosmic histories for four typical 
parameter cases of ‘no solution’ region. As seen below, the cosmic expansion rates evolve 
much differently depending on the parameter sets. 

(p, log|£o|) =(5, 65): As the time increases, the amplitudes of the third to fifth terms in 
eq. (2.13) rapidly increase with respect to that of the second term. The local minimal value 
in the function g(H;^,p) (eq. (2.13)) as a function of H then increases. At a certain time, 
the local minimal value becomes positive, and a real positive root for H disappears. We note 


10 





-6 -4 -2 0 2 4 6 8 10 12 14 


H (10“ 7 s -1 ) 

Figure 3. The functions g(H, £, p) as a function of H for the parameter set of (p, log |£o|> f(G i))=(5, 
65, 0). One line is the function at the cosmic time of t = 1.341 x 10 6 s which has a real positive 
solution marked by an open circle, while another is the function at t = 1.462 x 10 6 s after the real 
solution disappeared. 


that the amplitudes of the third to fifth terms are not larger than that of the second when the 
solution is lost although the ratios of amplitudes of former three terms and the latter term 
are rapidly increasing. 

Figure 3 shows the functions g(H,^,p) as a function of H for the parameter set of (p, 
log |£o|) f{G i))=(5, 65, 0). Two lines for t = 1.341 x 10 6 and t = 1.462 x 10 6 s are shown. The 
former time has a real positive solution marked by an open circle, while the latter corresponds 
to the time right after the real solution disappeared. 

Figure 4 shows respective terms in the function g(H , £, p) (eq. (2.13)) at t = 1.462 x 10 6 s 
as a function of H for the same parameter set as in figure 3. Dotted lines show the first and the 
second terms which exist in the standard cosmological model, while the solid (third term), 
dot-dashed (fourth), and dashed (fifth) lines show terms which exist only in the modified 
gravity model. This figure is for a time right after the proper solution is lost (figure 3). The 
vertical lines at H/( 1CD 7 s -1 ) ~ 7.6 correspond to asymptotic lines of —n 2 f{G) (fourth term) 
and k?£G (fifth). The loss of a solution is thus caused in the ‘no solution’ region by the three 
additional terms of the modified gravity in the function g(H,^,p). As seen in this case, a 
strong cancellation of the three terms can occur. 

(p, log |£o|) = (2, 84): In the epoch when the solution of H is lost in the BBN calculation, 
the second and fifth terms in eq. (2.13) are negative while the third and fourth terms are 
positive. The amplitudes of the modified-gravity terms satisfies |24J7 3 £| > \f(G)\ > |£G|. We 
note that relations in the standard model (eqs. (5.1)-(5.3)) are not satisfied in this case since 
the deviation in the expansion rate from standard case is large. Especially, the inequality 
G > 0 or -H < H 2 (eq. (2.6)) is satisfied in this epoch. The amplitude of the first term is 
negligible compared to other terms. At a certain time, the local minimal value in the function 
g(H;^,p) becomes larger than zero, and a real positive root of H disappears. 

Figure 5 shows the functions g(H,£,p) as a function of H for the parameter set of (p, 
l°g|fo|, /(Gi)) = (2, 84, 0). The functions are plotted for t = 1.193 x 10~ 2 s with a real 
positive solution marked by an open circle, and for t = 2.755 x 1CP 2 s after the real solution 
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Figure 4. Respective terms in the functions g(H,Gp) (eq. (2.13)) at t = 1.462 x 10 6 s as a function 
of H for the same parameter set as in figure 3. Dotted lines show the first and the second terms, 
while the solid, dot-dashed, and dashed lines show the third, fourth, and fifth terms, respectively. 



-0.8 -0.4 0 0.4 0.8 1.2 1.6 


H (10 2 s“ l ) 

Figure 5. The functions g(H, £, p) as a function of H for the parameter set of (p, log |£o|, /(Gj))=(2, 
84, 0). Two lines correspond to the functions at t = 1.193 x ICR 2 s with a real positive solution 
marked by an open circle and at t = 2.755 x 10“ 2 s after the real solution disappeared. 


disappeared. 

Figure 6 shows respective terms in the function g(H, £, p ) (eq. (2.13)) at t = 2.755 X ICR 2 
s as a function of H for the same parameter set as in figure 5. The lines correspond to the same 
quantities as in figure 4. The vertical lines at H/( 10 2 s” 1 ) ~ 0.6 correspond to asymptotic 
lines of — n 2 f(G) (fourth term) and k 2 £G (fifth). 

(p, log |^o|) =(—1.5, 78): When the solution of H is lost, the first and fourth terms are 
dominant, and are balanced with each other, and the second term, —6 H 2 , is much smaller 
than the two terms. The term k 2 /(G), however, becomes larger than the term 2 k 2 p with 
increasing time. The local maximal value of the g(H]^,p) function then becomes negative, 
and the proper root of H disappears. 
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Figure 6. Respective terms in the functions g(H,£,p) (eq. (2.13)) at t = 2.755 x 10~ 2 s as a function 
of H for the same parameter set as in figure 5. The lines correspond to the same quantities as in 
figure 4. 



Figure 7. The functions g{H 1 £,p) as a function of H for the parameter set of (p, log|£o|> 
/(Gi))=(— 1.5, 78, 0). Two lines correspond to the functions at t = 4.892 x 10~ 2 s with a real positive 
solution marked by an open circle and at t = 5.024 x 10 -2 s after the real solution disappeared. 


Figure 7 shows the functions g(H , £, p ) as a function of H for the parameter set of 
( p , log|£o|, 1.5, 78, 0). We plot the functions at t = 4.892 x 10 -2 s with a real 

positive solution marked by an open circle, and at t. = 5.024 X 10^ 2 s after the real solution 
disappeared. 

Figure 8 shows respective terms in the function g(H, £, p ) (eq. (2.13)) at t = 5.024 X 10~ 2 
s as a function of H for the same parameter set as in figure 7. The lines correspond to the 
same quantities as in figure 4. 

( p, log |£o|) =(—1-5, 82): In this case, a positive real solution exits. The solution, however, 
becomes much larger than that in SBBN, and this case is excluded based on the condition 
(4) (See section 3). 
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Figure 8. Respective terms in the functions g{H , £, p) (eq. (2.13)) at t = 5.024 x 10~ 2 s as a function 
of H for the same parameter set as in figure 7. The lines correspond to the same quantities as in 
figure 4. 



Figure 9. Cosmic expansion rates H as a function of Tg. The solid line corresponds to the present 
/(G) gravity model with (p, log |£ 0 |, f(G{))=(— 1.5, 82, 0), while the dashed line corresponds to SBBN. 
The difference between the two models becomes extremely large with decreasing temperature in the 
late epoch. 


Figure 9 shows the cosmic expansion rates H as a function of Tg. The solid line corre¬ 
sponds to the present f(G) gravity model with (p, log|£o| ; /(Gi))=(—1.5, 82, 0), while the 
dashed line corresponds to SBBN. When the expansion rate becomes 10 8 times larger than 
the rate in SBBN, the calculation is terminated. Although this case has a positive solution 
of H, this parameter region is safely excluded from the limits on elemental abundances. 

5.2.3 Effects on BBN 

Figure 10 shows the cosmic expansion rates H (s -1 ) as a function of Tg in the SBBN model 
(dashed line) and four cases of log|£o| = 74, 77, 80, and 83, respectively, (solid lines) with 
fixed parameters of p = 2 and f(G{) = 0 as an example case. We note that the third and fifth 
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Figure 10. Hubble expansion rate H as a function of Tg in SBBN (dashed line) and in the cases 
of log|£o| = 74, 77, 80, and 83, respectively, (solid lines). It is assumed that the power-law index is 
p = 2, and that the initial /(G) value is zero, i.e, /(G;) = 0. 



Figure 11. Mass fractions of H and 4 He (A' p and F p , respectively) and number ratios of other nuclides 
relative to H as a function of Tg. Solid lines show the abundances in the /(G) gravity model with the 
parameters ( p , log |£o|, [/(G)/(£G)]i)=(2, 74, 0), while dashed lines show those in SBBN. 


terms in eq. (2.13) linearly scale as the first standard term, i.e., —6 H 2 , in the case of p = 2 
as long as the relation H^l/t is not significantly broken (eqs. (5.6) and (5.8)). In SBBN, the 
relation of H ~ Tg is realized since the energy density is dominated by radiation in the BBN 
epoch. 

Standard BBN: Figure 11 shows evolutions of elemental abundances as a function of 
Tg in the SBBN model (dashed lines) and the f(G ) model with parameters ( p , log|£o|: 
[/(Cr)/(£Cr)]i) = (2, 74, 0) (solid lines). The quantities X p and Y p are mass fractions of 1 H and 
4 He, respectively. Other lines show the number ratios of nuclides to hydrogen, 1 H. 

Firstly, we review important reactions in SBBN model. At Tg > 10, by efficient weak 
reactions, the neutron-to-proton ratio is the equilibrium value of (n/p) eq = exp(—Q/T) with 
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Q = m n — m p = 1.293 MeV the mass difference. After the weak reaction freeze-out at 
Tg ~ 10, the /3-decay of the neutron is the dominant reaction for neutron destruction. When 
the temperature decreases to Tg < 1, the neutron is destroyed mainly via 1 H(n, 7 ) 2 H. D is 
predominantly produced via 4 H(n, 7 ) 2 H, and destroyed via 2 H(d, ra) 3 He and 2 H(d, p) 3 H. 3 H 
is produced via 2 H(d, p) 3 H and destroyed via 3 H(d, n) 4 He. 3 He is produced via 2 H(d, n) 3 He 
and destroyed via 3 He(n,p) 3 H. The primordial 3 H abundance is the sum of abundances of 3 H 
and 3 He produced during the BBN. Long after the BBN, 3 He nuclei /3-decay into 3 H nuclei 
with the half-life of 12.32 y. The final abundance of 3 He is larger than that of 3 H by about 
two orders of magnitude in SBBN. Therefore, the primordial 3 H abundance predominantly 
reflects the larger 3 He abundance during BBN. 

'Li is produced via 4 He(t, 7 ) 7 Li and destroyed via 7 Li(p, a) 4 He. 'Be is produced via 
4 He( 3 He, 7 ) 7 Be and destroyed via 'Be(n,p) 7 Li. The primordial 'Li abundance is the sum 
of abundances of 'Li and 'Be produced during the BBN. Long after the BBN, 'Be nuclei 
recombine with electrons and are transformed to ' Li nuclei via the electron capture process. 
'Be is produced more than 'Li in SBBN. The primordial 'Li abundance then predominantly 
reflects the larger 'Be abundance during BBN. (, Li is produced via 4 He(d, 7 ) 6 Li and destroyed 
via 6 Li(p, a) 3 He. 

Results of BBN in the f{G) model are compared with those in SBBN below. 

l°g |Co I = 74: The expansion rate becomes larger than that in SBBN at Tg < 0.4 (figure 10). 
The temperature then decreases in a time scale shorter than in SBBN after the temperature 
comes down to this critical value. Because of the faster expansion, the deuterium destruction 
is less efficient, and the D/H abundance is higher (figure 11). 

log I Co I = 77: Figure 12 shows elemental abundances as a function of Tg, similarly to figure 
11 , in the case of parameters ( p , log |Co|, [/(G)/(CG)]i)=(2, 77, 0) (solid lines). The expansion 
rate becomes larger than that in SBBN at Tg < 2 (figure 10). The cosmic time at a given 
temperature is then shorter. As a result, the n/p ratio is higher at the 4 He synthesis temper¬ 
ature of Tg ~ 1, and 4 He abundance is then slightly higher. We note that the 4 He synthesis 
occurs at slightly lower temperature than in SBBN because of the shorter expansion time 
scale for a given temperature. The destructions of D and 3 H freeze out earlier than in SBBN. 
The destructions are, therefore, less efficient, and the abundances D/H and 3 H/H are higher. 
The production and destruction rates of 3 He is higher because of higher D abundance. As a 
result, the abundance 3 He/H is slightly larger. Destructions of ‘Li and 6 Li freeze out earlier 
than in SBBN. The destructions are then less efficient, and the abundances 7 Li/H and 6 Li/H 
are higher. The 'Be abundance is very small because of efficient destruction by abundant 
neutrons via 'Be(n, p)‘ Li. 

l°g I Co I = 80: Figure 13 shows elemental abundances as a function of Tg, similarly to figure 
11, in the case of parameters ( p , log |£o|, [f (G) / (^G)]i) — (2, 80, 0) (solid lines). The expansion 
rate becomes larger than that in SBBN at Tg < 10. The 4 He production is incomplete, and 
the 4 He synthesis occurs at lower temperature than in SBBN. The 4 He abundance is then 
lower. The destructions of D and 3 H are less efficient, and the abundances D/H and 3 H/H are 
higher. The neutron abundance is very high since the expansion time scale is short compared 
to the neutron /3-decay half-life. 3 He is then effectively converted to 3 H by neutrons via 3 He(n, 
p) 3 H. The abundance 3 He/H is, therefore, smaller. The nuclei 7 Li and 6 Li are produced at 
lower temperature. The destructions are then less efficient, and the abundances 'Li/H and 
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Figure 12. Same as in figure 11, but solid lines correspond to the /(G) gravity model with the 
parameters (p, log|£ 0 |, [/(G)/(£G)]i)=(2, 77, 0). 



Figure 13. Same as in figure 11, but solid lines correspond to the /(G) gravity model with the 
parameters (p, log|£ 0 |, [/(G)/(£G)]i)=(2, 80, 0). 


6 Li/H are higher. The 'Be abundance is very small because of efficient destruction via 'Be(n, 
p) 7 Li. 

log |£o| = 83: Figure 14 shows elemental abundances as a function of Tg, similarly to figure 
11, in the case of parameters (p, log |£o|, [/(G)/(£G)]i)=(2, 83, 0) (solid lines). The expansion 
rate becomes larger than that in SBBN at Tg < 50. The 4 He production is incomplete, and 
the 4 He synthesis occurs at lower temperature than in SBBN. The 4 He abundance is then 
much lower. (Y p ~ 10~ 4 ). The destructions of D and 3 H are less efficient. Their productions 
are, however, also less efficient. Therefore, the abundances D/H and 3 H/H are higher than in 
SBBN although smaller than in case (d). The neutron abundance is very high. 3 He is then 
effectively converted to 3 H via 3 He(n, p) 3 H so that its abundance is smaller. The nuclei 'Li 
and 'Be are produced at lower temperature. Their abundances are, however, small because 
of less abundances of 4 He, 3 H, and 3 He. The 'Li destruction is less efficient. As a result, the 
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Figure 14. Same as in figure 11, but solid lines correspond to the /(G) gravity model with the 
parameters (p, log|£ 0 |, [/(G)/(£G)]i)=(2, 83, 0). 


abundance 'Li/H is smaller. The 6 Li abundance is very small because of a smaller abundance 
of 4 He. In addition, the 'Be abundance is effectively reduced by efficient destruction via 
7 Be(n, p) r Li. 


5.3 Case of £o > 0 and p = 2 
5.3.1 Parameter search 

Figure 15 shows the same boundary (black line) and contours of element abundances (colored 
lines) as in figure 2, but in the parameter plane of ([/(G)/(£G)]i, £o) for £o > 0- The power- 
law index is fixed as p = 2. When the condition [/(G)/(£G)]i < 2 is satisfied initially, 
the constraint of log£o 68 is derived from the requirement for the continuous existence 
of solution of H during the BBN. In the case of [/(G)/(£G)]i > 2, on the other hand, the 
constraint from the existence of solution disappears, and constraints from observational light 
element abundances are important. The reasons are described as follows. 

We suppose that the terms of modified gravity do not affect the cosmic expansion rate 
at the initial time very much. Then, because of G x < 0 (eqs. (2.6) and (5.3)) and £; > 0, 
the inequality ^G\ < 0 is satisfied. The region of /(G;)/(£iGi) < 0, therefore, corresponds to 
f{G\) > 0. As the value of /(G;) > 0 is larger, the negative-A-like effect of the f(G) term is 
stronger. As a result, smaller values of £o is allowed for smaller values of /(G;)/(£iGi) by the 
consideration of the existence of solution. On the other hand, the region of /(Gi)/(£;G;) > 0 
corresponds to /(G;) < 0. As the value of /(Gi) is smaller, it is more difficult to lose a proper 
solution of H. Therefore, the constraint from the existence of solution is weaker. 

Using physical quantities in standard model (as in eqs. (5.5)—(5.8) ), a somewhat detailed 
equation of the [f(G) — f{G\)] is derived as 


- k 2 [/(G) - /(G ; )] ~ 


(forp>4) 

Gi (for p < 4). 


(5.14) 


For example, when the power-law index is p = 2 as in the present case, it follows 


- At 2 [/(G) - /(Gi)] ~ 2« 2 6Gi < 0. 


(5.15) 
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Figure 15. Same boundary (black line) and contours of calculated light element abundances (colored 
lines) as in figure 2, but in the parameter plane of ([/(G)/(£G)];, £o) for £ 0 > 0. For this figure, the 
power-law index is fixed as p = 2. 


Then, when the inequality /(G;) < 2£iGi or /(Gi)/(£iGi) > 2 is satisfied, the condition 
—K?f(G) ~ k 2 [2^Gi — /(Gj)] > 0 always holds. The — k 2 /(G) term, therefore, operates 
similarly to a positive A. In this case, the real positive solution of H never disappears, and 
the ‘no solution’ region is absent in figure 15. In the region of [/(G)/(£G)]i > 2, when 
the amplitude £o is larger, the Hubble expansion rate is larger because of the /(G) term. 
Therefore, effects on the light element abundances are larger. Observational limits from 
elemental abundances then appear in the right upper corner in figure 15. For larger values of 
[/(G)/(£G)]i or smaller values of /(Gj), the increase of the cosmic expansion rate from that 
in SBBN is larger. Therefore, constraints from elemental abundances are severer. This fact 
explains the downward-sloping abundance contours in the region of [/(G)/(£G)]i > 2. 


p < 4: In general, a real positive solution of H exists when the following condition is satisfied, 


lim g(H-,£,p ) 

H ->+o v ' 


lim 2 k 2 p — 6 H 2 
H^+0 


-K 2 f(G)~ 12 k 4 (p + p) H 2 i 


lim 

H ->+0 


2 k 2 p - K 2 f(G ) - 6 [1 + 2k 4 (p + p) f] H 2 


2 k 2 P - k 2 /(G) > 0. 


(5.16) 


We note that in the limit as H approaches 0 from the right, the third term in eq. (2.13) 
(oc H a ) is negligible compared to the —6 H 2 term, and that the GB term has the limit value 
of limjy^+o G = —12 n 2 (p + p)H 2 (eqs. (2.6) and (2.14)). Equation (5.16) is satisfied when 
the following condition holds, 


/(G) ~ /(Gi) - —£iG; < 0 

4 — p 

/(GO > _J_ 

&G[ ~ 4-p’ 


(5.17) 


where eq. (5.14) was used in the first line. 
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Figure 16. Hubble expansion rate H as a function of Tg in SBBN (dashed line) and in the cases 
of log£o = 75, 76, and 78, respectively, with fixed parameters of p = 2 and [/(G)/(£G)]; = 4 (solid 
lines). 


p > 4: A real positive solution of H exists when eq. (5.16) is satisfied. In this case, however, 
eq. (5.14) is approximately true only when the cosmic expansion rate is not much different 
from that in the standard model. As the amplitude |/(G)| becomes relatively large with 
time, the expansion rate can be reduced. We define this epoch as the deviation time tdev 
Then, G oc H 2 (H + H 2 ) becomes very small after the deviation time. The contribution to 
the integration in the f{G ) function (eq. (3.1)) from the time after t dev is then small. As a 
result, the quantity of eq. (3.1) is mainly given by the value in the standard model at t = t dev . 
Equation (5.8), therefore, leads (cf. eq. (5.17)) to the condition, 


f(G) ~ /(<*) - ^4 ^dev < 0 

HGO y 4 
m dev ~p-4’ 


(5.18) 


where (£G)d ev is the value of £G at the deviation time. Before the deviation time, since the 
quantity (£G) has a simple scaling (eq. (5.8)), it follows that 


(£G)=£iGi (t/t if 4 . 

The following condition is then found, 

f(GQ ^ 4 /tdevV" 4 

(£iGi) — 4 V h ) 


(5.19) 


(5.20) 


5.3.2 Effects on BBN 

Figure 16 shows the cosmic expansion rates H (s _1 ) as a function of Tg in the SBBN model 
(dashed line) and three cases of log£o = 75, 76, and 78, respectively, (solid lines) with fixed 
parameters of p = 2 and [/(G)/(£G)]; = 4. 
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Figure 17. Same as in figure 11, but solid lines correspond to the /(G) gravity model with the 
parameters (p, log£ 0 , [/(G)/(£G)]i)=(2, 75, 4). 


(log£o> [/(G)/(£G)]i) = (75, 4): Figure 17 shows elemental abundances as a function of Tg, 
similarly to figure 11, in the case of parameters (p, log£o, [/(G)/(£Cr)]i)=(2, 75, 4) (solid lines). 
The expansion rate becomes larger than that in SBBN at Tg < 0.7 (figure 16). Deuterium 
destruction is then less efficient, and the final abundance D/H is higher. 3 H abundance is also 
somewhat larger. The neutron abundance remains to be high since the temperature decreases 
effectively without the change of cosmic time after the expansion rate deviates from that in 
SBBN. Destructions of 'Li and (, Li, and 'Be production freeze out earlier than in SBBN. 
Therefore, the abundances 'Li/H and h Li/H are higher while the 'Be/H abundance is lower. 

(log£o> [/(G)/(£G)];) = (76, 4): Figure 18 shows elemental abundances as a function of Tg, 
similarly to figure 11, in the case of parameters (p, log£o, [/(G)/(£G)];) = (2, 76, 4) (solid 
lines). The expansion rate becomes larger than that in SBBN at Tg < 1 (figure 16). In 
this case, destructions of D, 3 H, and 3 He are less efficient, and the final abundances of D/H, 
3 H/H, and 3 He/H are higher than those in SBBN. Also, destructions of 'Li and 6 Li, and 
'Be production freeze out earlier than in SBBN, similarly to the case (a). It then results in 
significantly higher 'Li/H and () Li/H and lower ‘Be/H values. 

(log^o? [/(G)/(£G)]i) = (78, 4); Figure 19 shows elemental abundances as a function of Tg, 
similarly to figure 11, in the case of parameters ( p , log£o> [/(G)/(£G)]i) = (2, 78, 4) (solid 
lines). The expansion rate becomes larger than that in SBBN at Tg < 3. The cosmic time at 
a given temperature is then shorter. As a result, the n/p ratio is higher at the 4 He synthesis 
temperature of Tg ~ 1, and 4 He abundance is then slightly higher. We note that the 4 He 
synthesis occurs at slightly lower temperature than in SBBN because of the shorter expansion 
time scale for a given temperature. The destructions of D and 3 H are less efficient, and the 
abundances D/H and 3 H/H are higher. The neutron abundance remains to be high (n/H 
~ 10~ 2 ). 3 He is then effectively converted to 3 H by abundant neutrons via 3 He(n, p) 3 H. The 
abundance 3 He/H is, therefore, smaller. The nuclei 'Li and 6 Li are produced at slightly lower 
temperature than in SBBN. Since the 'Li destruction is less efficient at lower temperature, the 
abundance 'Li/H is higher. The 'Be abundance is very small because of efficient destruction 
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Figure 18. Same as in figure 11, but solid lines correspond to the /(G) gravity model with the 
parameters (p, log£ 0 , [/(G)/(£G)]i)=(2, 76, 4). 



Figure 19. Same as in figure 11, but solid lines correspond to the /(G) gravity model with the 
parameters (p, log£ 0 , [/(G)/(£G)]i)=(2, 78, 4). 


by neutrons via 'Be(n, p)'Li. The abundance l, Li/H is higher than in SBBN. It is, however, 
reduced at T 9 ~ 0.5 — 0.1 via 6 Li(n, a) 3 H, which is never important in SBBN. 

5.4 Case of £0 < 0 and p = 2 
5.4.1 Parameter search 

Figure 20 shows the same contours as in figure 15 but for £0 < 0. Parameter regions in figures 
15 and 20 looks mirror symmetric with respect to the f (G {)/= 2 line. The solution 
boundary and abundance contours are explained as follows. 

We again suppose that the terms of modified gravity do not affect the cosmic expansion 
rate at the initial time very much. The inequalities G\ < 0 and £; < 0 are satisfied, and the 
inequality £;Gi > 0 holds. The region of f (G\) / (£,\G{) > 0, therefore, corresponds to f(G\) > 
0. Larger values of f(G{) then lead to stronger effects of the negative-A-like f(G) term. 
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Figure 20. Same as in figure 15, but in the parameter plane of ([/(G)/(£G)]i, |Co|) for Co < 0. 


Smaller values of Co are then allowed for larger values of /(G;)/(CiGi) by the consideration 
of the existence of solution. On the other hand, the region of /(G;)/(CiGi) < 0 corresponds 
to f(Gi) < 0. When f(G\) is smaller, it is more difficult to lose a proper solution of H, and 
the constraint from the existence of solution is weaker. However, the cosmic expansion rate 
is increased more, and light element abundances are affected more. As a result, constraints 
from the observational limits on the abundances are stronger, which can be seen from the 
upward-sloping abundance contours for [f(G)/(£G )]; < 2 . 

Similarly to the case of Co > 0, eqs. (5.14) and (5.19) are true before the deviation time. 


p < 4: The condition for the existence of proper solution of H is given by f(G ) < 0. This 
is satisfied if 


m) ^ 4 
CiGi ~4 -p 


(5.21) 


p > 4: The condition for the existence of the solution H is satisfied when the following 
condition holds, 


f(Gj) ^ 4 (t dev \ p ~ 4 

(CiGi) ~p-4 V ti ) 


(5.22) 


For example, for the power-law index p = 2, the condition for a proper solution is given 
(eq. (5.21)) by /(Gj)/(CiGj) < 2. If this condition is not met, i.e., /(Gj)/(CiG;) > 2, the term 
f(G) > 0 induces a negative-A effect, and the solution of H can be lost. 


5.4.2 Effects on BBN 

Figure 21 shows the cosmic expansion rates H (s _1 ) as a function of T 9 in the SBBN model 
(dashed line) and three cases of log |Co I = 75, 76, and 78, respectively, with fixed parameters 
of p = 2 and [/(G)/(CG)]i = 4 (solid lines). 

Trends of nucleosynthesis are very similar to those in section 5.3.2. 

(log |CoI? [/(G)/(CG)];) = (75, 4): Figure 22 shows elemental abundances as a function of Tg, 
similarly to figure 11 , in the case of parameters (p, log|Co|j [/(G)/(CG)];) = ( 2 , 75, 4) (solid 
lines). The expansion rate becomes larger than that in SBBN at Tg < 0.7 (figure 21). For 
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Figure 21. Hubble expansion rate H as a function of Tg in SBBN (dashed line) and in the cases 
of log |£ 0 | — 75, 76, and 78, respectively, with fixed parameters of p = 2 and [/(G)/(£G)]i = 4 (solid 
lines). 



Figure 22. Same as in figure 11, but solid lines correspond to the /(G) gravity model with the 
parameters (p, log|£ 0 |, [/(G)/(£G)]i)=(2, 75, 4). 


the same reasons as in case (a) in section 5.3.2, the abundances n/H, D/H, 3 He/H, 'Li/H, 
and 6 Li/H are higher while 'Be/H is lower than in SBBN. 

(log |£o|, [/(G)/(£Cr)]i) = (76, 4): Figure 23 shows elemental abundances as a function of Tg, 
similarly to figure 11, in the case of parameters (p, log|£o|, [/(G)/(£G)]j)=(2, 76, 4) (solid 
lines). The expansion rate becomes larger than that in SBBN at Tg < 1 (figure 21). Although 
the trends of abundances are qualitatively the same as those in the case (a), differences in 
abundances from the SBBN are larger in case (b) because of larger difference in the expansion 
rate. 

(log |£o|, [/(G)/(£G)]i) = (78, 4): Figure 24 shows elemental abundances as a function of Tg, 
similarly to figure 11, in the case of parameters (p, log|£o|, [/(G)/(£G)]i)=(2, 78, 4) (solid 
lines). The expansion rate becomes larger than that in SBBN at Tg < 3 (figure 21). Because 
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Figure 23. Same as in figure 11, but solid lines correspond to the /(G) gravity model with the 
parameters (p, log|£ 0 |, [/(G)/(£G)]i)=(2, 76, 4). 



Figure 24. Same as in figure 11, but solid lines correspond to the /(G) gravity model with the 
parameters (p, log|£ 0 |, [/(G)/(£G)]i)=(2, 78, 4). 


of the faster expansion, the n/p ratio is higher at Tg ~ 1, and 4 He abundance is slightly 
higher. The 4 He synthesis occurs at slightly lower temperature than in SBBN. For the same 
reasons as in case (c) in section 5.3.2, the abundances n/ H, D/H, 'Li/H, and 6 Li/H are higher 
while 3 He/H, and "Be/H are lower than in SBBN. 

6 Summary 

We studied effects of the /(G) gravity on BBN. It was assumed that a /(G) term exists in the 
BBN epoch of Tg = [100, 0.01]. The functional form was taken to be £ = df /dG = £o(i/io) p 
where to = 1 s is a typical time scale, and £o and p are parameters. Under this assumption, 
the /(G) model can be described with three parameters, i.e., the coefficient £o> the power-law 
index p, and the initial value of /(G;). We then showed a method to solve physical variables 
during BBN consistently taking account of the modified cosmic expansion rate. 
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Numerical calculations of BBN in the /(G) model were performed, and results were ana¬ 
lyzed. It was then found that a proper solution for the cosmic expansion rate during the BBN 
epoch does not exist in some parameter region. In addition, we compared calculated results 
of primordial light element abundances with observational data. Constraints on parameters 
of the f(G ) gravity are then derived from the observed abundances. Since a change in the 
cosmic expansion rate easily deviates primordial D abundance from observational limits, the 
observed D abundance gives the strongest constraint on the parameters (figures 2, 15, 20). 
We thus obtained allowed parameter regions where the universe expands properly and final 
abundances of light elements in the BBN epoch are consistent with observations. 

In the case of 6) > 0 and /(G;) = 0 (section 5.1), the existence of the solution for the 
cosmic expansion rate determines the allowed parameter region of 6) and p predominantly. 
The constraint is derived as 


(71-^ (3.5 - p) (for p < 3.5) 

log6)< { (71-70) (for 3.5 <p< 4.2) (6.1) 

l 70 - 6 (p-4.2) (for 4.2 < p). 

In the case of 6) < 0 and f(G\) = 0 (section 5.2), on the other hand, the parameter region 
is constrained from the limit on elemental abundances for p < 3.8 and from the existence of 
the solution for 3.8 < p. The derived constraint is 

f 74- |(3.8 -p) (for p < 3.8) 

log |£o| < { 74 - 12.5 (p - 3.8) (for 3.8 < p < 4) (6.2) 

[ 71.5 - 6.5 (p- 4) (for 4 < p). 


In the case of 6) > 0 and p = 2 (section 5.3), we derived a constraint on the parameter region 
of 6) and [/(G)/(£G)]i which is the ratio of the quantities /(G ) and (£G) at the initial time of 
computation. The constraint comes from the existence of the solution for [/(G)/(£G)]i < 2, 
and from the limit on elemental abundances for 2 < [/(G)/(£G)];. The derived constraint is 


log 6) ^ 


68 (for [/(G)/(£G)]i < 2) 
71 (for 2 < [/(G)/(£G)]i). 


(6.3) 


In the case of 6) < 0 and p = 2 (section 5.4), in contrast, the parameter region is constrained 
from the limit on elemental abundances for [/(G)/(£G)]i < 2 and from the existence of the 
solution for 2 < [/(G)/(£G)] ;. The derived constraint is 


log |6)| < 


71 (for [/(G)/(£G)]i < 2) 
68 (for 2 < [/(G)/(£G)]i). 


(6.4) 


These are constraints on the /(G) term during the BBN epoch. 
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